{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pulse-based Variational Quantum Eigensolver Algorithm\n",
    "\n",
    "\n",
    "*Copyright (c) 2021 Institute for Quantum Computing, Baidu Inc. All Rights Reserved.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Outline\n",
    "**Note: Running the program of this tutorial may cost more than 50 credit points in your Quantum Hub account.**\n",
    "\n",
    "This tutorial introduces how to implement the variational quantum eigensolver algorithm at the pulse level. The outline of this tutorial is as follows:\n",
    "- Introduction to variational quantum eigensolver\n",
    "- Introduction to Pulse-Based variational quantum eigensolver\n",
    "- Preparation\n",
    "- Construct Hamiltonian\n",
    "- Optimize two-qubit gates\n",
    "- Construct problem Hamiltonian\n",
    "- Pulse-based ansatz and optimization\n",
    "- Summary"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to variational quantum eigensolver\n",
    "\n",
    "Variational Quantum Eigensolver (VQE) is a popular quantum algorithm for approximating the ground state energy of molecules on Noisy Intermediate-Scale Quantum (NISQ) computers. VQE is mainly used to estimate the minimum eigenvalue of a given Hamiltonian and find its ground state. For near-term quantum computers, critical gate errors, decoherence, and poor connectivity limit the depth of quantum circuits. However, the VQE algorithm can work with low-depth quantum circuits. Hence, the VQE algorithm is considered an ideal candidate to utilize NISQ devices to solve real-world problems.\n",
    "\n",
    "In VQE, the essential task is to calculate the ground state energy of a given discretized molecular Hamiltonian $\\hat{H}_{\\rm mole}$ by a parametrized trail state $|\\psi(\\vec{\\theta})\\rangle$. The trail state $|\\psi(\\vec{\\theta})\\rangle$ is generated by a parameterized ansatz. At the same time, we employ the classical optimization methods to find a set of optimal $\\vec{\\theta}$ to minimize the expectation value $E = \\langle \\psi(\\vec{\\theta}) | \\hat{H}_{\\rm mole} | \\psi(\\vec{\\theta}) \\rangle$, which represents the approximate ground state energy $E_0$ of the Hamiltonian $\\hat{H}_{\\rm mole}$:\n",
    "$$\n",
    "E_0 = {\\rm min}_{\\vec{\\theta}} \\langle \\psi(\\vec{\\theta}) | \\hat{H}_{\\rm mole} | \\psi(\\vec{\\theta}) \\rangle\n",
    "$$\n",
    "\n",
    "In this tutorial, we will introduce the method to implement VQE for solving the ground state energy of a hydrogen molecule $H_2$. We simulate the VQE on a superconducting system at the pulse level while considering multiple nonideal factors. The following figure shows a common ansatz for generating trail states:\n",
    "\n",
    "![VQE](figures/vqe-circuit.png)\n",
    "\n",
    "In this example, we no longer use CNOT gates, because it is not directly implemented in the superconducting platform. However, we use the superconducting hardware-efficient two-qubit gates, i.e., cross-resonance (CR) gates, which also can realize maximal entanglements using single-qubit gates. Thus, the above four-qubit circuit is formed by single-qubit gates and CR gates. The matrix of an ideal CR gate is:\n",
    "\n",
    "$$\n",
    "\n",
    "\\hat{U}_{\\rm CR}(\\alpha) = \\begin{bmatrix}\n",
    "\\cos{\\frac{\\alpha}{2}} & -i\\sin{\\frac{\\alpha}{2}} & 0 & 0 \\\\\n",
    "-i\\sin{\\frac{\\alpha}{2}} & \\cos{\\frac{\\alpha}{2}} & 0 & 0 \\\\ \n",
    "0 & 0 & \\cos{\\frac{\\alpha}{2}} & i\\sin{\\frac{\\alpha}{2}} \\\\\n",
    "0 & 0 & i\\sin{\\frac{\\alpha}{2}} & \\cos{\\frac{\\alpha}{2}} \n",
    "\\end{bmatrix}.\n",
    "\n",
    "$$\n",
    "\n",
    "Here, we choose $\\alpha = -\\frac{\\pi}{2}$. More details about the CR gate can be found [here](https://quanlse.baidu.com/#/doc/tutorial-cr). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction to Pulse-Based variational quantum eigensolver\n",
    "\n",
    "In this tutorial, we study the VQE algorithm at the pulse level, which we call **pulse-based VQE**. Unlike the classical VQE which optimizes the rotation parameters of the single-qubit gates from the logical quantum circuit, pulse-based VQE directly takes the pulse parameters as the optimization parameters to find the minimal loss value (i.e. ground state energy). The following figure shows the difference between the pulse-based VQE and the standard VQE:\n",
    "\n",
    "![VQE](figures/vqe-scheme.png)\n",
    "\n",
    "To implement pulse-based VQE, we need to translate the logical quantum circuit into a **pulse-based quantum circuit**, i.e., the logical rotating gates $R_x(\\theta)$ and $R_y(\\theta)$ are replaced by the corresponding control pulses on $X$ and $Y$ channels respectively with variable amplitudes - which we call **pulse-based gates**:\n",
    "\n",
    "![VQE](figures/vqe-translate.png)\n",
    "\n",
    "In this figure, $U_{\\rm ENT}$ is a unitary matrix which can create the required entanglement - the details will be explained in the following sections. Here, we use a new notation representing the parameters of **pulse-based gates**:\n",
    "$$\n",
    "\\vec{A} = [A_1, \\cdots, A_m, \\cdots, A_M],\n",
    "$$\n",
    "where $M$ is the total number of **pulse-based gates**; $A_m$ denotes the amplitude of the Gaussian waveform of the $m$-th **pulse-based gate**:\n",
    "$$\n",
    "\\Omega_m(t) = A_m e^{-(\\frac{t - \\tau_m}{\\sqrt{2} \\sigma_m}) ^2}.\n",
    "$$\n",
    "Other pulse parameters such as the width $\\sigma_m$ and the central position $\\tau_m$ will be fixed throughout the procedure. Based on the introduction of **pulse-based VQE**, we no longer need to generate an optimized drive pulse according to the parameters from the logical quantum circuit. This improves the efficiency of VQE and the accuracy of the results significantly.\n",
    "\n",
    "In the above sections, we give a brief introduction to the standard and pulse-based VQE. In the following sections, we will demonstrate the implementation of pulse-based VQE using Quanlse step by step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preparation\n",
    "After you have successfully installed Quanlse, you could run the Quanlse program below following this tutorial. To run this particular tutorial, you would need to import the following packages from Quanlse and other commonly-used Python libraries:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "from numpy import linalg, min, random, trace, dot, savez, load, identity, kron\n",
    "\n",
    "from math import pi\n",
    "from functools import reduce\n",
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# This module creates the Hamiltonian dictionary\n",
    "from Quanlse.Utils import Hamiltonian as qham\n",
    "\n",
    "# These functions help us perform matrix calculations\n",
    "from Quanlse.Utils.Tools import unitaryInfidelity, tensor\n",
    "\n",
    "# These functions define useful operator matrices\n",
    "from Quanlse.Utils.Operator import sigmaX, sigmaY, sigmaZ, sigmaI\n",
    "\n",
    "# This function generates wave data\n",
    "from Quanlse.Utils.Waveforms import makeWaveData\n",
    "\n",
    "# This function uploads jobs to Quanlse Cloud Service and receives results\n",
    "from Quanlse.remoteSimulator import remoteSimulatorRunHamiltonian as runHamiltonian\n",
    "\n",
    "# This module defines matrices of the frequently used quantum gates\n",
    "from Quanlse.QOperation import FixedGate\n",
    "\n",
    "# This module saves the PBVQE results\n",
    "from Quanlse.Define import outputPath\n",
    "\n",
    "# Generate the path of npz file\n",
    "localFile = os.path.join(outputPath, f'pbvqe.npz')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct Hamiltonian\n",
    "\n",
    "First, we define some necessary constants, including the sampling period of the arbitrary wave generator (AWG), the number of qubits in the circuit, and the system's energy levels to consider."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sampling period (Nano second)\n",
    "dt = 2.0\n",
    "\n",
    "# Number of qubits\n",
    "qubits = 4\n",
    "\n",
    "# System energy level\n",
    "level = 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we define the hardware parameters of the superconducting qubits. The items in the `freq` list are $\\omega_{\\rm q0}, \\omega_{\\rm q1}, \\omega_{\\rm q2}, \\omega_{\\rm q3}$; the items in the `coupling` list save the coupling information of qubit 0-1, 1-2, 2-3, 3-0, respectively. By using the rotating wave approximation, we define the system in a rotaing frame with the qubits' frequencies at: $\\omega_{\\rm RWA} = \\omega_{\\rm q0} = \\omega_{\\rm q2} = 4.914 \\times 2\\pi$ GHz."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the hardware parameters of the qubits (GHz)\n",
    "freq = [4.914 * (2 * pi), 5.114 * (2 * pi), 4.914 * (2 * pi), 5.114 * (2 * pi)]\n",
    "\n",
    "# Define the coupling strength (GHz)\n",
    "coupling = [\n",
    "    [[0, 1], 0.016 * (2 * pi)],\n",
    "    [[1, 2], 0.016 * (2 * pi)],\n",
    "    [[2, 3], 0.016 * (2 * pi)],\n",
    "    [[3, 0], 0.016 * (2 * pi)]\n",
    "]\n",
    "\n",
    "# Frequency of rotating frame (GHz)\n",
    "rwa = 4.914 * (2 * pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we fix the gate time for all of the single-qubit and two-qubit gates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Gate time (Nano second)\n",
    "tg_2q = 200\n",
    "tg_1q = 64"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we create the Hamiltonian according to the following hardware structure: each qubit couples with its neighbors and the coupling strength is always constant.\n",
    "\n",
    "![VQE](figures/vqe-topo_structure.png)\n",
    "\n",
    "Here, we construct the Hamiltonian of the system above:\n",
    "$$\n",
    "\\hat{H}_{\\rm total} = \\sum_{q=0}^{3} \\delta_{q} \\hat{a}^{\\dagger}_{q}\\hat{a}_{q} + \\frac{1}{2}\\sum_{q=0}^{3}(\\hat{a}_{q}\\hat{a}^{\\dagger}_{q+1}+\\hat{a}^{\\dagger}_{q}\\hat{a}_{q+1}) + \\sum_{q=0}^{3}\\Omega_{q}^x (t) \\hat{\\sigma}_{q}^{x} + \\sum_{q=0}^{3}\\Omega_{q}^y (t) \\hat{\\sigma}_{q}^{y} + \\sum_{q=0}^{3}\\Omega_{q}^z (t) \\hat{\\sigma}_{q}^{z} ,\n",
    "$$\n",
    "where $\\hat{a}_{q}$ and $\\hat{a}^{\\dagger}_{q}$ are the annihilation and creation operators on the $q$-th qubit respectively. $\\hat{\\sigma}^x_{q}, \\hat{\\sigma}^y_{q}$ and $\\hat{\\sigma}^z_{q}$ are the Pauli operators on the $q$-th qubit. $\\delta_{q}=\\omega_{q} - \\omega_{\\rm RWA}$ is the detuning of the $q$-th qubit, $\\Omega_q^{x,y,z}(t)$ are the envelope functions of the drive pulses or of the magnetic flux control. In Quanlse, the Hamiltonian above can be constructed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the Hamiltonian\n",
    "vqeham = qham.createHam(title=\"3q-2l-vqe\", dt=dt, qubitNum=qubits, sysLevel=level)\n",
    "\n",
    "# Add the coupling terms\n",
    "for item in coupling:\n",
    "    q0, q1 = item[0][0], item[0][1]\n",
    "    qham.addCoupling(vqeham, f\"coupling{q0}-{q1}\", onQubits=[q0, q1], g=item[1] / 2)\n",
    "\n",
    "for qubit in range(qubits):\n",
    "    # Add the detuning terms\n",
    "    detuning = freq[qubit] - rwa\n",
    "    qham.addDrift(vqeham, name=f\"detuning{qubit}\", onQubits=qubit, matrices=sigmaZ(), amp=detuning)\n",
    "    # Add the control terms\n",
    "    qham.addControl(vqeham, name=f\"x{qubit}\", onQubits=qubit, matrices=sigmaX())\n",
    "    qham.addControl(vqeham, name=f\"y{qubit}\", onQubits=qubit, matrices=sigmaY())\n",
    "    qham.addControl(vqeham, name=f\"setdetuning{qubit}\", onQubits=qubit, matrices=sigmaZ())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For more information on constructing Hamiltonian using Quanlse, please see [tutorial-single-qubit-gate](https://quanlse.baidu.com/#/doc/tutorial-single-qubit)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Optimize two-qubit gates\n",
    "\n",
    "In this example, we use cross-resonance gates as the entangling gates. For more information on the CR gate, please see [tutorial-cr-gate](https://quanlse.baidu.com/#/doc/tutorial-cr). Due to the direct coupling of the neighboring qubits, applying X pulses on one qubit will end up affecting two qubits. Hence, we need to take this factor into account when we design our pulses in order to suppress the unwanted effects. \n",
    "\n",
    "![VQE](figures/vqe-crosstalk.png)\n",
    "\n",
    "Here, we use `qham.subSystem()` to extract two three-qubit sub-systems from `vqeham` for optimizing the CR gates - one of which is qubits 0-1-2 and the other is qubits 1-2-3. On these sub-systems, we set $\\hat{U}_{\\rm goal}=I \\otimes \\hat{U}_{\\rm CR}$ as the goal unitary matrix to optimize the pulses' parameters - the goal is to generate a CR gate on the second and the third qubits of the sub-system.\n",
    "\n",
    "We define a function `make_cr_pulse()` to define the waveforms. The Gaussian microwave drive pulses are applied to the second qubits of the sub-systems. We fix the width and the center position of the Gaussian wave while setting the amplitude as the first parameter to optimize. The second parameter to be optimized is the amplitude of the flux detuning (square wave) applied on the first qubit. Note that the other `setdetuning` terms are used to transform the rotating frame to a particular frequency."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_cr_pulse(ham, sub_sys_3q, drive_freq, amp, shift, t0, t):\n",
    "    \"\"\" Assemble the pulses for CR gates \"\"\"\n",
    "    wave_list = []\n",
    "    cr_para = {\"a\": amp, \"tau\": tg_2q / 2, \"sigma\": tg_2q / 8}\n",
    "    wave_list.append(makeWaveData(ham, f\"x{sub_sys_3q[1]}\", f=\"gaussian\", t0=t0, t=t, para=cr_para))\n",
    "    # frame transformation\n",
    "    wave_list.append(makeWaveData(ham, f\"setdetuning{sub_sys_3q[0]}\", f=\"square\", t0=t0, t=t, para={\"a\": rwa - drive_freq + shift}))\n",
    "    wave_list.append(makeWaveData(ham, f\"setdetuning{sub_sys_3q[1]}\", f=\"square\", t0=t0, t=t, para={\"a\": rwa - drive_freq}))\n",
    "    wave_list.append(makeWaveData(ham, f\"setdetuning{sub_sys_3q[2]}\", f=\"square\", t0=t0, t=t, para={\"a\": rwa - drive_freq}))\n",
    "    return wave_list"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we define a function `optimize_cr` to perform the optimization. The optimal parameters will be saved for further usage."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimize_cr(sub_sys_3q, drive_freq):\n",
    "    \"\"\" Realize a CR gate on the second & third qubits \"\"\"\n",
    "\n",
    "    crham = qham.subSystem(vqeham, sub_sys_3q, \"3q-2l-cropt\")\n",
    "    ugoal = tensor([identity(2), FixedGate.CR.getMatrix()])\n",
    "\n",
    "    def cr_loss(x):\n",
    "        # Clear and add waves\n",
    "        qham.clearWaves(crham)\n",
    "        # Generate and add waves for CR gate implementation\n",
    "        wave_list = make_cr_pulse(crham, sub_sys_3q, drive_freq, x[0], x[1], 0, tg_2q)\n",
    "        qham.addWaveData(crham, wave_list)\n",
    "        # Simulate the system's evolution and obtain the infidelity\n",
    "        unitary = qham.getUnitary(crham)\n",
    "        infidelity = unitaryInfidelity(ugoal, unitary, 3)\n",
    "        return infidelity\n",
    "\n",
    "    opt = optimize.dual_annealing(cr_loss, [(-2, 2), (-0.2, 0.2)], maxiter=60)\n",
    "    print(\"Min infidelity:\", opt[\"fun\"])\n",
    "    return opt[\"x\"][0], opt[\"x\"][1]\n",
    "\n",
    "lhl_1q_x, lhl_0q_z = optimize_cr([0, 1, 2], 4.914 * 2 * pi)\n",
    "hlh_1q_x, hlh_0q_z = optimize_cr([1, 2, 3], 5.114 * 2 * pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Construct problem Hamiltonian\n",
    "\n",
    "In this section, we introduce the method for estimating the ground state energy for hydrogen molecule $H_2$ at the pulse level. Here we will skip the details on fermion-to-qubit mapping, please visit [Paddle Quantum](https://github.com/PaddlePaddle/Quantum/blob/master/tutorial/VQE/VQE_EN.ipynb) for more information. First, we define a function `pauli_str_to_matrix()` to convert the Pauli string to the Hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def pauli_str_to_matrix(pauli_str, n):\n",
    "    \"\"\"\n",
    "    Convert the Pauli string in Hamiltonian\n",
    "    \"\"\"\n",
    "    def NKron(AMatrix, BMatrix, *args):\n",
    "        return reduce(\n",
    "            lambda result, index: kron(result, index),\n",
    "            args,\n",
    "            kron(AMatrix, BMatrix), )\n",
    "    pauli_dict = {'i': sigmaI(), 'x': sigmaX(), 'y': sigmaY(), 'z': sigmaZ()}\n",
    "    # Parse pauli_str; 'x0,z1,y4' to 'xziiy'\n",
    "    new_pauli_str = []\n",
    "    for coeff, op_str in pauli_str:\n",
    "        init = list('i' * n)\n",
    "        op_list = op_str.split(',')\n",
    "        for op in op_list:\n",
    "            pos = int(op[1:])\n",
    "            assert pos < n, 'n is too small'\n",
    "            init[pos] = op[0]\n",
    "        new_pauli_str.append([coeff, ''.join(init)])\n",
    "\n",
    "    # Convert new_pauli_str to matrix; 'xziiy' to NKron(x, z, i, i, y)\n",
    "    matrices = []\n",
    "    for coeff, op_str in new_pauli_str:\n",
    "        sub_matrices = []\n",
    "        for op in op_str:\n",
    "            sub_matrices.append(pauli_dict[op])\n",
    "        if len(op_str) == 1:\n",
    "            matrices.append(coeff * sub_matrices[0])\n",
    "        else:\n",
    "            matrices.append(coeff * NKron(sub_matrices[0], sub_matrices[1], *sub_matrices[2:]))\n",
    "\n",
    "    return sum(matrices)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use the Hamiltonian of a $H_2$ molecule with an interatomic distance of $d=74$ pm. The following data is generated from [Paddle Quantum](https://github.com/PaddlePaddle/Quantum/blob/master/tutorial/VQE/VQE_EN.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "target_ham = [\n",
    "    [-0.042078976477822, 'i0'],\n",
    "    [ 0.177712874651399, 'z0'],\n",
    "    [ 0.177712874651399, 'z1'],\n",
    "    [-0.242742805131446, 'z2'],\n",
    "    [-0.242742805131462, 'z3'],\n",
    "    [ 0.170597383288005, 'z0,z1'],\n",
    "    [ 0.044750144015351, 'y0,x1,x2,y3'],\n",
    "    [-0.044750144015351, 'y0,y1,x2,x3'],\n",
    "    [-0.044750144015351, 'x0,x1,y2,y3'],\n",
    "    [ 0.044750144015351, 'x0,y1,y2,x3'],\n",
    "    [ 0.122933050561837, 'z0,z2'],\n",
    "    [ 0.167683194577189, 'z0,z3'],\n",
    "    [ 0.167683194577189, 'z1,z2'],\n",
    "    [ 0.122933050561837, 'z1,z3'],\n",
    "    [ 0.176276408043195, 'z2,z3']\n",
    "]\n",
    "H_matrix = pauli_str_to_matrix(target_ham, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The theoretical ground state energy is:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the theoretical eigenvalue\n",
    "eig_val, eig_state = linalg.eig(H_matrix)\n",
    "min_eig_H = min(eig_val.real)\n",
    "print(f\"Ground state energy: {min_eig_H} Ha\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Pulse-based ansatz and optimization\n",
    "\n",
    "We design the pulse-based VQE ansatz using one of the most commonly-used templates for a $H_2$ molecule. The following figure shows one layer of the ansatz, containing 3 single-qubit gates for each qubit. Each single-qubit gate holds one parameter as the amplitude of a Gaussian envelope, in which the width and the center position are fixed.\n",
    "\n",
    "![VQE](figures/vqe-scheduling.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we define a function `makeWaveSchedule()` to schedule the pulse sequence according to the pulse-based ansatz shown above. Argument `x` is a list of optimization parameters (i.e., pulse parameters $\\vec{A}$); `job_waves` is a list of waveform data generated by `makeWaveData()`, which saves the detailed information of the user-defined waveform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def makeWaveSchedule(x):\n",
    "    \"\"\" Generate waves for pulse-based circuit \"\"\"\n",
    "    job_waves = []\n",
    "    qham.clearWaves(vqeham)\n",
    "    start_time = 0\n",
    "    depth = int(len(x) / 12)\n",
    "    for d in range(depth):\n",
    "        start_time = d * (4 * tg_2q + 3 * tg_1q)\n",
    "        # Add pulses for single-qubit gates\n",
    "        for g in range(4):\n",
    "            # X control\n",
    "            pulse_para_x = {\"a\": x[12 * d + g], \"tau\": tg_1q / 2, \"sigma\": tg_1q / 8}\n",
    "            job_waves.append(makeWaveData(vqeham, f\"x{g}\", t0=start_time, t=tg_1q, f=\"gaussian\", para=pulse_para_x))\n",
    "            # Y control\n",
    "            pulse_para_y = {\"a\": x[12 * d + 4 + g], \"tau\": tg_1q / 2, \"sigma\": tg_1q / 8}\n",
    "            job_waves.append(makeWaveData(vqeham, f\"y{g}\", t0=start_time + tg_1q, t=tg_1q, f=\"gaussian\", para=pulse_para_y))\n",
    "            # X control\n",
    "            pulse_para_x = {\"a\": x[12 * d + 8 + g], \"tau\": tg_1q / 2, \"sigma\": tg_1q / 8}\n",
    "            job_waves.append(makeWaveData(vqeham, f\"x{g}\", t0=start_time + 2 * tg_1q, t=tg_1q, f=\"gaussian\", para=pulse_para_x))\n",
    "            # Set detuning\n",
    "            job_waves.append(makeWaveData(vqeham, f\"setdetuning{g}\", t0=start_time, t=tg_1q * 3, f=\"square\", para={\"a\": rwa - freq[g]}))\n",
    "        \n",
    "        # Add pulses for two-qubit gates\n",
    "        cr_start_time = start_time + tg_1q * 3\n",
    "        wave_list = make_cr_pulse(vqeham, [3, 0, 1], 5.114 * 2 * pi, hlh_1q_x, hlh_0q_z, cr_start_time, t=tg_2q)\n",
    "        job_waves.extend(wave_list)\n",
    "        wave_list = make_cr_pulse(vqeham, [0, 1, 2], 4.914 * 2 * pi, lhl_1q_x, lhl_0q_z, cr_start_time + tg_2q * 1, t=tg_2q)\n",
    "        job_waves.extend(wave_list)\n",
    "        wave_list = make_cr_pulse(vqeham, [1, 2, 3], 5.114 * 2 * pi, hlh_1q_x, hlh_0q_z, cr_start_time + tg_2q * 2, t=tg_2q)\n",
    "        job_waves.extend(wave_list)\n",
    "        wave_list = make_cr_pulse(vqeham, [2, 3, 0], 4.914 * 2 * pi, lhl_1q_x, lhl_0q_z, cr_start_time + tg_2q * 3, t=tg_2q)\n",
    "        job_waves.extend(wave_list)\n",
    "    \n",
    "    return job_waves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we use the gradient-based optimization method (L-BFGS-B) provided by `Scipy`. In this method, information regarding the gradient for all parameters must be provided. We use two-point finite difference method to approximate the gradient:\n",
    "$$\n",
    "\\frac{\\partial{\\rm Loss}(\\vec{A})}{\\partial A_m} = \\frac{{\\rm Loss}(A_0, \\cdots, A_m + \\epsilon, \\cdots, A_{M-1}) - {\\rm Loss}(A_0, \\cdots, A_m - \\epsilon, \\cdots, A_{M-1})}{2\\epsilon} ,\n",
    "$$\n",
    "where $\\vec{A} = [A_0, \\cdots, A_{M-1}]$ is the pulse parameter list, $\\epsilon$ is a small number, and the ${\\rm Loss}(\\vec{A})$ is defined as,\n",
    "$$\n",
    "{\\rm Loss}(\\vec{A}) =  \\langle \\psi(\\vec{A}) | \\hat{H}_{\\rm mole} | \\psi(\\vec{A}) \\rangle.\n",
    "$$\n",
    "Trail state $\\psi(\\vec{A})$ is generated by the pulse-based ansatz. Finite difference requires a large number of samples. For example, given the number of pulse parameter $M$, we need $2M$ samples to estimate the approximated gradient. Hence, we use Quanlse Cloud Service to accelerate this procedure. \n",
    "\n",
    "To utilize Quanlse Cloud Service, we need to import `Define` and acquire a token, which can be obtained on [Baidu Quantum-Hub](http://quantum-hub.baidu.com)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the loss function\n",
    "import copy\n",
    "from Quanlse import Define\n",
    "Define.hubToken = \"\"\n",
    "loss_history = []"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then, we define the VQE `loss` function. In this function, we simulate the evolution of the pulse-based circuit at $\\vec{x}$ (the amplitude parameters of pulse-based gates) and compute the gradient at this point by the finite difference method we mentioned above. In each iteration, we input the current pulse parameter list $\\vec{x}$ into the loss function and package all samples into `waveList`. `waveList` contains $2M$ samples for solving the gradient and $1$ sample for obtaining the loss value. \n",
    "\n",
    "Finally, we integrate all the waves into one list, i.e., `waveList`, and submit it to Quanlse Cloud Service by function `runHamiltonian()`. After about 15 to 20 seconds, we receive the result, which will be saved to the `Output` folder as a JSON file. At the same time, the variable `result` will be assigned a list that contains all the simulation results corresponding to the `waveList`.\n",
    "\n",
    "**Note:** Each item of `waveList` contains all the waves for pulse-based VQE, generated by the function `makeWaveSchedule()` we have defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def loss(x):\n",
    "    global loss_history\n",
    "    # Add wave for current point\n",
    "    waveList = []\n",
    "    waveList.append(makeWaveSchedule(x))\n",
    "    \n",
    "    # Add wave for calculating gradient\n",
    "    for xid in range(len(x)):\n",
    "        xlist = copy.deepcopy(x)\n",
    "        xlist[xid] -= 1e-8\n",
    "        waveList.append(makeWaveSchedule(xlist))\n",
    "        xlist[xid] += 2 * 1e-8\n",
    "        waveList.append(makeWaveSchedule(xlist))\n",
    "\n",
    "    # Simulate the evolution\n",
    "    result = runHamiltonian(vqeham, jobList=waveList)\n",
    "\n",
    "    # Calculate the loss function\n",
    "    loss_list = []\n",
    "    for item in result:\n",
    "        unitary = item[\"unitary\"]\n",
    "        loss_val = (unitary.conj().T @ H_matrix @ unitary).real[0][0]\n",
    "        loss_list.append(loss_val)\n",
    "    \n",
    "    # Calculate the gradients\n",
    "    gradient = []\n",
    "    for index in range(len(x)):\n",
    "        gradient.append((loss_list[2 + index * 2] - loss_list[1 + index * 2]) / (1e-8 * 2))\n",
    "    \n",
    "    print(\"Loss function:\", loss_list[0])\n",
    "    loss_history.append(loss_list[0])\n",
    "    return loss_list[0], gradient"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we use `fmin_l_bfgs_b()` provided by `Scipy` to minimize the loss function, which corresponds to the ground state energy. We start with randomly initialized parameters. In this example, we set the ansatz depth $d=3$, define the bounds for all pulse parameters $\\vec{x}$ and limit the maximal numbers of iteration.\n",
    "\n",
    "**Note**: this optimization may take longer than 15 minutes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "depth = 3\n",
    "\n",
    "loss_history = []\n",
    "\n",
    "init_paras = [random.rand() for _ in range(depth * 12)]\n",
    "bounds = [(-1.5, 1.5) for _ in range(depth * 12)]\n",
    "x, f, d = optimize.fmin_l_bfgs_b(loss, init_paras, fprime=None, bounds=bounds, maxiter=100)\n",
    "\n",
    "# Save the loss history to a file for further usage\n",
    "savez(localFile, loss_history)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"The estimated ground state energy is: {f} Ha\")\n",
    "print(\"Total iteration:\", d[\"nit\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we see, the optimization converges at high precision, and the final number of iteration is 77.\n",
    "Then, we plot the optimized energy values to the iteration number:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load the loss_history list from the npz file.\n",
    "loss_history = load(localFile)['arr_0']\n",
    "\n",
    "# Plot the figures\n",
    "import matplotlib.pyplot as plt\n",
    "plt.plot(range(len(loss_history)), loss_history, label=\"Energy\")\n",
    "plt.axhline(min_eig_H, c=\"gray\", ls=\"--\", lw=1.0)\n",
    "plt.xlabel(\"Iteration\")\n",
    "plt.ylabel(\"Energy (Ha)\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can plot the pulse sequence using the `plotWaves()` function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Print the waveforms.\n",
    "qham.addWaveData(vqeham, makeWaveSchedule(x))\n",
    "qham.plotWaves(vqeham, color=['red', 'green', 'blue'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary\n",
    "After reading this tutorial on pulse-based VQE, the users are encouraged to explore other advanced research which are different from this tutorial.\n",
    "\n",
    "## References\n",
    "\n",
    "\\[1\\] [Peruzzo, Alberto, et al. \"A variational eigenvalue solver on a photonic quantum processor.\" *Nature communications* 5 (2014): 4213.](https://doi.org/10.1038/ncomms5213)\n",
    "\n",
    "\\[2\\] [Moll, Nikolaj, et al. \"Quantum optimization using variational algorithms on near-term quantum devices.\" *Quantum Science and Technology* 3.3 (2018): 030503.](https://doi.org/10.1088/2058-9565/aab822)\n",
    "\n",
    "\\[3\\] [Kandala, Abhinav, et al. \"Hardware-efficient variational quantum eigensolver for small molecules and quantum magnets.\" *Nature* 549.7671 (2017): 242-246.](https://doi.org/10.1038/nature23879)\n",
    "\n",
    "\\[4\\] [Rigetti, Chad, and Michel Devoret. \"Fully microwave-tunable universal gates in superconducting qubits with linear couplings and fixed transition frequencies.\" *Physical Review B* 81.13 (2010): 134507.](https://doi.org/10.1103/PhysRevB.81.134507)\n",
    "\n",
    "\\[5\\] [Meitei, Oinam Romesh, et al. \"Gate-free state preparation for fast variational quantum eigensolver simulations: ctrl-VQE.\" *arXiv preprint arXiv:2008.04302 (2020)*.](https://arxiv.org/abs/2008.04302)\n",
    "\n",
    "\\[6\\] [Wilhelm, Frank K., et al. \"An introduction into optimal control for quantum technologies.\" *arXiv preprint arXiv:2003.10132* (2020).](https://arxiv.org/abs/2003.10132)\n",
    "\n",
    "\\[7\\] [Krantz, Philip, et al. \"A quantum engineer's guide to superconducting qubits.\" *Applied Physics Reviews* 6.2 (2019): 021318.](https://aip.scitation.org/doi/abs/10.1063/1.5089550)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
